Analysis for Dowry Deaths

load needed libraries:

library(dplyr)        # A package for data manipulation
library(sf)           # Simple feature for R
library(spdep)        # Functions and tests for evaluating spatial patterns 
library(tidyr)        # Tools to create tidy data
library(INLA)         # Integrated Nested Laplace Approximation package
library(ggplot2)      # A package for creating maps and graphs
library(viridis)      # A package providing color palettes 
library(patchwork)    # A package to compose plots

# For tables in RMarkdown
library(knitr)
library(kableExtra)

1. Import and explore the data

Import the file with the data and call the data.frame object as data.

load("~/Desktop/AA project/DS1_CrimeUttarPradesh/CrimeUttarPradesh.RData")

Compute the total number of cases of dowry death per year, and format the output in a table. To display a table with additional formatting, we can use the knitr::kable function:

kable(data %>%
        group_by(year) %>%
         summarise(dowry_observed = sum(dowry), dowry_expected=sum(e_dowry)), booktabs = T, caption = "dowry death cases in Uttar Pradesh by year") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center")

dowry death cases in Uttar Pradesh by year
year	dowry_observed	dowry_expected
2001	2210	1708.608
2002	1893	1760.031
2003	1322	1811.454
2004	1708	1862.877
2005	1563	1914.300
2006	1798	1965.723
2007	2074	2017.146
2008	2237	2068.568
2009	2231	2119.992
2010	2201	2171.414
2011	2321	2222.837
2012	2242	2274.260
2013	2332	2325.683
2014	2468	2377.106

Import the shape file carto_up using the function st_read from sf package and call the object as carto_up_sf.

carto_up_sf <- st_as_sf(carto_up)

Then, plot the spatial object carto_up_sf using ggplot2 package.

ggplot() + 
  geom_sf(data = carto_up_sf, color = "lightblue", fill = "white") +
  coord_sf() +    #axis limits and CRS
  theme_bw() +    # dark-on-light theme
  theme(axis.title = element_text(size = 14),
        axis.text = element_text(size = 12))

2. Spatial model

In order to use the same data structure for both the space-only model and later the space-time model, a new set of data is formed by aggregating both the observed and expected counts over time. A Poisson-log linear model is then fitted, assuming a BYM2 model for the random effects. Let each areal unit $i$ be indexed by the integers $1, 2,\dots,N$. The model is specified as:

\[ \begin{eqnarray} O_i &\sim & \text{Poisson}(\rho_i E_i)\\ \eta_i &= & \log \rho_i = b_0 + b_i\\ \boldsymbol{b} &= & \frac{1}{\sqrt{\tau_b}}(\sqrt{1-\phi}\boldsymbol{v}_{*} + \sqrt{\phi}\boldsymbol{u}_{*})\\ \end{eqnarray} \]

where $\boldsymbol{v}_{*}$ and $\boldsymbol{u}_{*}$ are standardised versions of $\boldsymbol{u}$ and $\boldsymbol{v}$. In detail, we specify the spatial random effect $\boldsymbol{b}$ using a re-parameterisation of the Besag-York-Molli'e prior [@Besag1991], which is a convolution of an intrinsic conditional autoregressive (CAR) model and an independent and identically distributed Gaussian model, where

$u_i$ is the spatially structured component defined by an intrinsic CAR prior: $\boldsymbol{u}\sim ICAR(\boldsymbol{W}, \sigma^2_u)$,
$v_i$ the unstructured component defined with prior: $v_s \overset{iid}{\sim} \text{Normal}(0,\sigma^2_v)$.

Define the neighbors and create the weights list.

carto_up_sf_nb = poly2nb(carto_up_sf, queen=TRUE)
summary(carto_up_sf_nb)

## Neighbour list object:
## Number of regions: 70 
## Number of nonzero links: 324 
## Percentage nonzero weights: 6.612245 
## Average number of links: 4.628571 
## 2 regions with no links:
## 4 11
## Link number distribution:
## 
##  0  1  2  3  4  5  6  7  8  9 
##  2  2  2 13 13 16 12  6  3  1 
## 2 least connected regions:
## 45 60 with 1 link
## 1 most connected region:
## 16 with 9 links

Convert the list of neighbors to inla format using the function nb2WB(). call the graph as carto_india_sf.adj

nb2INLA("carto_up_sf.graph", carto_up_sf_nb)
carto_up_sf.adj=paste(getwd(), "/carto_up_sf.graph", sep ="")

Aggregate observed and expected dowry deaths cases over the geographical areas ID_area.

data_agg = data %>% group_by (ID_area) %>%
  summarize(dowry_observed = sum(dowry),
            dowry_expected = sum(e_dowry)) %>% 
  dplyr::rename(dowry_obs = dowry_observed, dowry_exp = dowry_expected)

compute the standardized morbidity ratios for the dowry cases (SMRs)

data_agg = data_agg %>% mutate (dowry_SMR = dowry_obs/dowry_exp)

Produce a spatial map of the aggregated SMRs using ggplot2 package. For the map use the following breakpoints [min,0.4], (0.4-0.6], (0.6-0.8], (0.8,1], (1,1.2], (1.2-1.4], (1.4-1.6], (1.6-max].

Remember that, before to produce the map, you need to join the sf object and the data frame data_agg. To do so you can use the function left_join from the library dplyr.

data_agg$dowry_SMRcat = cut(data_agg$dowry_SMR, 
                      breaks=c(min(data_agg$dowry_SMR), 
                               0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 
                               max(data_agg$dowry_SMR)), include.lowest = T)

map_SMR = left_join(carto_up_sf, data_agg, by = c("ID_area" = "ID_area"))

and plot:

ggplot() + geom_sf(data = map_SMR, col = NA) + aes(fill = dowry_SMRcat) +
  theme_bw() + scale_fill_viridis_d() + 
  guides(fill=guide_legend(title="SMR for dowry deaths cases"))

Map of the average dowry deaths SMRs over the period 2001-2014 in Uttar Pradesh

Fit the hierarchical Poisson log-linear model in INLA. Remember to create an ID for the areas, i.e. index vector, and use a bym2 prior for the spatial random effect setting the PC priors [@Simpson2017] seen in Session 2.2. Monitor also the Watanabe-Akaike information criterion (WAIC) including control.compute=list(waic=TRUE)

ID = seq(1,70)
formula_BYM2 = dowry_obs ~ f(ID, model="bym2", graph="carto_up_sf.graph",
        hyper=list(prec = list(
        prior = "pc.prec",
        param = c(0.5 / 0.31, 0.01)),
        phi = list(
        prior = "pc",
        param = c(0.5, 2 / 3))))  

sBYM.model = inla(formula=formula_BYM2, family="poisson", data=data_agg, E=E, control.compute=list(waic=TRUE))

## Warning in set.warn("E", nm): Argument 'E=E' expanded to NULL or gave an error.
##   This might be an error and you are requested to check this out.
##   Move on with default values...

## Warning in inla.model.properties.generic(inla.trim.family(model), mm[names(mm) == : Model 'bym2' in section 'latent' is marked as 'experimental'; changes may appear at any time.
##   Use this model with extra care!!! Further warnings are disabled.

Obtain the posterior summary statistics (mean and posterior probability that the residual is above 1 - (or log-residual is above 0)) of the parameters of interest:

#Relative risks
RR_sBYM = c()

for(i in 1:70){
  RR_sBYM[i] = inla.emarginal(function(x) exp(x), 
        sBYM.model$marginals.random$ID[[i]])
}

#Posterior probabilities
RR_sBYM_marg = sBYM.model$marginals.random$ID[1:70]
PP_sBYM = lapply(RR_sBYM_marg, function(x) {1-inla.pmarginal(0,x)})

Obtain the posterior estimates from the spatial model to be plotted, that is (i) the area level posterior mean of the residual RRs and (ii) the posterior probability (PP) that the residual RRs > 1.

resRR_PP = data.frame(resRR=RR_sBYM, 
                       PP=unlist(PP_sBYM),
                      ID_area=data_agg[,1])

Using ggplot2 package, produce a map of the posterior mean of the residual RRs and the posterior probabilities that the residual RRs are > 1

For the map of the posterior mean of the residual RRs, use the following breakpoints [min,0.4], (0.4-0.6], (0.6-0.8], (0.8,1], (1,1.2], (1.2-1.4], (1.4-1.6], (1.6-max].

resRR_PP$resRRcat = cut(resRR_PP$resRR, breaks=c(min(resRR_PP$resRR), 
                  0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 
                  max(resRR_PP$resRR)),include.lowest = T)

For the map of the probabilities PP that the residual RRs is > 1, use the following breakpoints [0-0.2], (0.2-0.8], (0.8-1].

# breakpoints
resRR_PP$PPcat = cut(resRR_PP$PP, c(0, 0.2, 0.8, 1.00), include.lowest = TRUE)

Remember to join sf object and data frame with the posterior estimates

map_RR_PP = left_join(carto_up_sf, resRR_PP, by = c("ID_area" = "ID_area"))

Produce the maps of the posterior mean of the residual RRs and the posterior probabilities PP using ggplot2 package.

p1 = ggplot() + geom_sf(data = map_RR_PP) + aes(fill = resRRcat) +
  theme_bw() + scale_fill_brewer(palette = "PuOr") + 
  guides(fill=guide_legend(title="RR")) + ggtitle("RR Spatial model for dowry death cases") + 
  theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))

p2 = ggplot() + geom_sf(data = map_RR_PP) + aes(fill = PPcat) +
  theme_bw() +
  scale_fill_viridis(
    option = "plasma", name="PP",
    discrete = T,
    direction = -1,
    guide = guide_legend(
      title.position = 'top',
      reverse = T
    )) +  ggtitle("PP Spatial model for dowry death") + theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")) 

p1|p2

3. Estimate the spatial fraction

As the BYM2 has the structured (CAR) and unstructured (iid) components it might be useful to get some ideas about the strength of the spatially structured components as this would indicate the level of clustering in the data. To do so we ca simply obtain the posterior summary of the phi hyperparameter. How do you interpret it?

sBYM.model$summary.hyperpar

##                       mean        sd 0.025quant  0.5quant 0.975quant      mode
## Precision for ID 3.6652700 0.6959124  2.4619005 3.6101026  5.1934901 3.5114492
## Phi for ID       0.4734963 0.1716023  0.1642017 0.4682986  0.8071051 0.4445385

ANSWER: This tells us that about 1/2 of the spatial variability is explained by the spatially structured component - which makes sense if we look at the map of the RR, which shows a degree of spatial clustering.

4. Spatio-temporal model (no interaction)

Now, we extend the above analysis to a separable space-time model without interactions. For the temporal component, we use the specification introduces in the lecture with a temporal unstructured random effect and a structured one (RW1 prior).

Let each areal unit $i$ be indexed by the integers $1, 2,\dots,N$. As in the spatial case, we use a Poisson distribution to model the number of hospital admission $O_it$, in area $i$ at time $t$. The mathematical specification of the model includes now an additional temporal dependence term, which can be modeled using a non-stationary random walk prior: ${i} ({i-1}, ^2_{}).
The model implement in this practical assumes no space-time interaction and a spatial convolution with random walk in time:

\[ \begin{eqnarray} O_{it} &\sim & \text{Poisson}(\rho_{it} E_{it}) \\ \log \rho_{it} &= & b_0 + b_i + \gamma_t + \psi_t \\ \boldsymbol{b} &= & \frac{1}{\sqrt{\tau_b}}(\sqrt{1-\phi}\boldsymbol{v}_{*} + \sqrt{\phi}\boldsymbol{u}_{*})\\ \gamma_t & \sim & \hbox{RW(1)}\\ \psi_t & \sim & N(0,\sigma^2_{\psi})\\ \end{eqnarray} \]

As seen before $\boldsymbol{v}_{*}$ and $\boldsymbol{u}_{*}$ are standardised versions of $\boldsymbol{u}$ and $\boldsymbol{v}$. In detail, we specify the spatial random effect $\boldsymbol{b}$ using a re-parameterisation of the Besag-York-Molli'e prior, which is a convolution of an intrinsic conditional autoregressive (CAR) model and an independent and identically distributed Gaussian model, where

$u_i$ is the spatially structured component defined by an intrinsic CAR prior: $\boldsymbol{u}\sim ICAR(\boldsymbol{W}, \sigma^2_u)$,
$v_i$ the unstructured component defined with prior: $v_s \overset{iid}{\sim} \text{Normal}(0,\sigma^2_v)$.

First prepare the data, joining in the shapefile to make sure that the order is the same and then create an ID for time and one for space.

#Join the data with the shapefile so the order of the shapefile is maintained.  
data_ST = left_join(carto_up_sf, data, by="ID_area")

#Rename the columns of Observed and Expected as we did before
data_ST = data_ST  %>% dplyr::rename(dowry_obs = dowry, dowry_exp = e_dowry)

#Create the ID for year (time)
data_ST$ID.time = data_ST$year - 2001

#Create the ID for space
data_ST$ID.space = rep(seq(1,70),each=14)

Run the model in INLA, monitoring the WAIC and using PC priors (as specified in Session 2.2) for the bym2 and rw1 models (note: for the rw1 model you could try the following specification f(ID.time,model="rw1", hyper=list(prec = list(prior = "pc.prec",param = c(0.5 / 0.31, 0.01))))). Call the output as stBYM.model

formula_ST_noint = dowry_obs ~ f(ID.space, model="bym2", graph=carto_up_sf.adj,
                            hyper=list(prec = list(
                            prior = "pc.prec",
                            param = c(0.5 / 0.31, 0.01)),
                            phi = list(
                            prior = "pc",
                            param = c(0.5, 2 / 3)))) + f(ID.time, model="rw1", hyper=list(prec = list(
                            prior = "pc.prec",
                            param = c(0.5 / 0.31, 0.01))))
                            
stBYM.model = inla(formula=formula_ST_noint, family="poisson", data=data_ST, E=E, 
                     control.compute=list(waic=TRUE))

## Warning in set.warn("E", nm): Argument 'E=E' expanded to NULL or gave an error.
##   This might be an error and you are requested to check this out.
##   Move on with default values...

Create the posterior mean for the spatial (as we did in point 10) and temporal effects

#Spatial Relative risks
RR_stBYM = c()

for(i in 1:70){
  RR_stBYM[i] = inla.emarginal(function(x) exp(x), 
        stBYM.model$marginals.random$ID.space[[i]])
}

#Posterior probabilities (for spatial RR)
RR_stBYM_marg = stBYM.model$marginals.random$ID.space[1:70]
PP_stBYM = lapply(RR_stBYM_marg, function(x) {1-inla.pmarginal(0,x)})

#Temporal Relative risks and CI95
RR_stRW_RR = c()
RR_stRW_lo = c()
RR_stRW_hi = c()

for(i in 1:14){
  #Posterior mean
  RR_stRW_RR[i] = inla.emarginal(function(x) exp(x), 
        stBYM.model$marginals.random$ID.time[[i]])
  #2.5% quantile 
  RR_stRW_lo[i] = inla.qmarginal(0.025,inla.tmarginal(function(x) exp(x), stBYM.model$marginals.random$ID.time[[i]]))
  #97.5% quantile 
  RR_stRW_hi[i] = inla.qmarginal(0.975, inla.tmarginal(function(x) exp(x), stBYM.model$marginals.random$ID.time[[i]]))
}

RR_stRW = data.frame(RR=RR_stRW_RR,low=RR_stRW_lo,high=RR_stRW_hi)

Plot the temporal residual RRs (RR_stWR)

Temp1 = ggplot(RR_stRW, aes(seq(2001,2014), RR)) + geom_line() + ggtitle("ST model No Int") + geom_ribbon(aes(ymin=low,ymax=high), alpha=0.2) + labs(x="year")

Temp1

Map the spatial residual RRs (RR_stBYM) with ggplot2 package using the following breakpoints [min,0.4], (0.4-0.6], (0.6-0.8], (0.8,1], (1,1.2], (1.2-1.4], (1.4-1.6], (1.6-max]. Compare this map against the map of the residual RR obtained from the spatial model.

resRR_PP_st = data.frame(resRR=RR_stBYM, 
                       PP=unlist(PP_stBYM),
                      ID_area=data_agg[,1])
# breakpoints
resRR_PP_st$resRRcat = cut(resRR_PP_st$resRR, breaks=c(min(resRR_PP_st$resRR), 
                  0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 
                  max(resRR_PP_st$resRR)),include.lowest = T)

resRR_PP_st$PPcat = cut(resRR_PP_st$PP, c(0, 0.2, 0.8, 1.00), include.lowest = TRUE)

map_RR_ST = left_join(carto_up_sf, resRR_PP_st, by = c("ID_area" = "ID_area"))

p3 = ggplot() + geom_sf(data = map_RR_ST) + aes(fill = resRRcat) +
  theme_bw() + scale_fill_brewer(palette = "PuOr") + 
  guides(fill=guide_legend(title="RR")) +  ggtitle("RR ST model") +
  theme(text = element_text(size=15), 
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")) 

p4 = ggplot() + geom_sf(data = map_RR_ST) + aes(fill = PPcat) +
  theme_bw() +
  scale_fill_viridis(
    option = "plasma",
    name = "PP ST model",
    discrete = T,
    direction = -1,
    guide = guide_legend(
      title.position = 'top',
      reverse = T
    )) +  ggtitle("PP ST model") + theme(text = element_text(size=15), 
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))

(p1|p2) / (p3|p4)

Spatio-temporal model: Map of the residual RRs and posterior probabilities for dowry death cases

5. Spatio-temporal model (type I interaction)

Now, we extend the above spatio-temporal analysis to a separable space-time model with type I interaction:

\[ \begin{eqnarray} O_{it} &\sim & \text{Poisson}(\rho_{it} E_{it}) \\ \log \rho_{it} &= & b_0 + b_i + \gamma_t + \psi_t + \delta_{it} \\ \boldsymbol{b} &= & \frac{1}{\sqrt{\tau_b}}(\sqrt{1-\phi}\boldsymbol{v}_{*} + \sqrt{\phi}\boldsymbol{u}_{*})\\ \gamma_t & \sim & \hbox{RW(1)}\\ \psi_t & \sim & N(0,\sigma^2_{\psi})\\ \delta_{it} & \sim & \hbox{Normal}(0, \sigma^2_{\delta}) \end{eqnarray} \]

Specify the formula and run the model in INLA. Remember that you need to create an index which goes from 1 to the length of the dataset (i.e. the space x time). For the iid model defining the interaction term, use the PC prior previously used for the rw1 model. Call the output stIntI.BYM.model (remember to monitor the WAIC)

data_ST$ID.space.time = seq(1,dim(data_ST)[1])

formula_ST_intI = dowry_obs ~ f(ID.space, model="bym2", graph=carto_up_sf.adj,
                            hyper=list(prec = list(
                            prior = "pc.prec",
                            param = c(0.5 / 0.31, 0.01)),
                            phi = list(
                            prior = "pc",
                            param = c(0.5, 2 / 3)))) + 
                      f(ID.time,model="rw1", hyper=list(prec = list(
                            prior = "pc.prec",
                            param = c(0.5 / 0.31, 0.01))))+
                      f(ID.space.time,model="iid", hyper=list(prec = list(
                            prior = "pc.prec",
                            param = c(0.5 / 0.31, 0.01))))
                    
                            
stIntI.BYM.model = inla(formula=formula_ST_intI, family="poisson", data=data_ST, E=E,
                        control.compute=list(dic=TRUE, waic=TRUE))

## Warning in set.warn("E", nm): Argument 'E=E' expanded to NULL or gave an error.
##   This might be an error and you are requested to check this out.
##   Move on with default values...

Create the posterior mean for the spatial and temporal effects and compare with the ST model results without interaction

#Spatial Relative risks
RR_stIntI.BYM = c()

for(i in 1:70){
  RR_stIntI.BYM[i] = inla.emarginal(function(x) exp(x), 
        stIntI.BYM.model$marginals.random$ID.space[[i]])
}

#Posterior probabilities (for spatial RR)
RR_stIntI.BYM_marg = stIntI.BYM.model$marginals.random$ID.space[1:70]
PP_stIntI.BYM = lapply(RR_stIntI.BYM_marg, function(x) {1-inla.pmarginal(0,x)})

#Temporal Relative risks and CI95
RR_stIntI.RW_RR = c()
RR_stIntI.RW_lo = c()
RR_stIntI.RW_hi = c()

for(i in 1:14){
  #Posterior mean
  RR_stIntI.RW_RR[i] = inla.emarginal(function(x) exp(x), 
        stIntI.BYM.model$marginals.random$ID.time[[i]])
  #2.5% quantile 
  RR_stIntI.RW_lo[i] = inla.qmarginal(0.025,inla.tmarginal(function(x) exp(x), stIntI.BYM.model$marginals.random$ID.time[[i]]))
  #97.5% quantile 
  RR_stIntI.RW_hi[i] = inla.qmarginal(0.975, inla.tmarginal(function(x) exp(x), stIntI.BYM.model$marginals.random$ID.time[[i]]))
}

RR_stIntI.RW = data.frame(RR=RR_stIntI.RW_RR,low=RR_stIntI.RW_lo,high=RR_stIntI.RW_hi)

Plot the temporal residual RRs (RR_stWR)

Temp2 = ggplot(RR_stIntI.RW, aes(seq(2001,2014), RR)) + geom_line() + ggtitle("ST model Int I") + geom_ribbon(aes(ymin=low,ymax=high), alpha=0.2) + labs(x="year")

Temp1 | Temp2

Map the spatial residual RRs (RR_stIntI.BYM) with ggplot2 package using the following breakpoints [min,0.4], (0.4-0.6], (0.6-0.8], (0.8,1], (1,1.2], (1.2-1.4], (1.4-1.6], (1.6-max]. Compare this map against the map of the residual RR obtained from the spatio-temporal model with no interaction. Try to comment the results

resRR_PP_stIntI = data.frame(resRR=RR_stIntI.BYM, 
                       PP=unlist(PP_stIntI.BYM),
                      ID_area=data_agg[,1])
# breakpoints
resRR_PP_stIntI$resRRcat = cut(resRR_PP_stIntI$resRR, breaks=c(min(resRR_PP_stIntI$resRR), 
                  0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 
                  max(resRR_PP_stIntI$resRR)),include.lowest = T)

resRR_PP_stIntI$PPcat = cut(resRR_PP_stIntI$PP, c(0, 0.2, 0.8, 1.00), include.lowest = TRUE)

map_RR_ST.IntI = left_join(carto_up_sf, resRR_PP_stIntI, by = c("ID_area" = "ID_area"))

p5 = ggplot() + geom_sf(data = map_RR_ST.IntI) + aes(fill = resRRcat) +
  theme_bw() + scale_fill_brewer(palette = "PuOr") + 
  guides(fill=guide_legend(title="RR")) +  ggtitle("RR ST model Int I") +
  theme(text = element_text(size=15), 
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")) 

p6 = ggplot() + geom_sf(data = map_RR_ST.IntI) + aes(fill = PPcat) +
  theme_bw() +
  scale_fill_viridis(
    option = "plasma",
    name = "PP ST model Int I",
    discrete = T,
    direction = -1,
    guide = guide_legend(
      title.position = 'top',
      reverse = T
    )) +  ggtitle("PP ST model Int I") + theme(text = element_text(size=15), 
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))

(p1|p2) / (p3|p4) / (p5|p6)

Spatio-temporal model: Map of the residual RRs and posterior probabilities for dowry death cases

ANSWER: We basically see that across the different models there is no difference in the spatial residuals. Let’s now look at the space-time interaction.

Plot the space-time interaction

data_ST$intI = stIntI.BYM.model$summary.random$ID.space.time$mean

data_ST$intI_cat = cut(data_ST$intI,  breaks=c(-1,-0.05, 
                  -0.01, 0.01, 0.05, 1),include.lowest = T)
ggplot() +
  geom_sf(data = data_ST, aes(fill = intI_cat))+ theme_bw() +  scale_fill_brewer(palette = "PuOr") + 
            guides(fill=guide_legend(title=NULL)) + 
            theme(text = element_text(size=20), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank()) +
facet_wrap(~ year, ncol = 3, labeller=labeller(ID.year=c("1"="2001","2"="2002","3"="2003","4"="2004","5"="2005", "6"="2006", "7"="2007", "8"="2008", "9"="2009", "10"="2010", "11"="2011", "12"="2012", "13"="2013", "14"="2014"))) +
labs("")

There is difference but I do not see a clear pattern

Get a table of the hyperparameters. How do you interpret this table?

dat.hyper2 = 
  round(
  data.frame(median = stIntI.BYM.model$summary.hyperpar[,4],
    LL = stIntI.BYM.model$summary.hyperpar[,3], 
    UL = stIntI.BYM.model$summary.hyperpar[,5]),
  digits = 3)

row.names(dat.hyper2) = 
  rownames(stIntI.BYM.model$summary.hyperpar)

knitr::kable(dat.hyper2, caption = "Posterior median and 95% CrI of hyperparameters for dowry death cases.") %>%  kable_styling(bootstrap_options = "striped", full_width = F, position = "center")

Posterior median and 95% CrI of hyperparameters for dowry death cases.
	median	LL	UL
Precision for ID.space	3.641	2.495	5.230
Phi for ID.space	0.471	0.170	0.811
Precision for ID.time	48.913	19.384	113.376
Precision for ID.space.time	37.605	30.397	46.542

Compare the WAIC from the three different models: spatial, spatio-temporal without interaction and spatio-temporal with interaction. What do you observe?

dat.WAIC = data.frame(model = c("Spatial", "SpatTemp no int", "SpatTemp typeI"), 
                       WAIC = round(c(sBYM.model$waic$waic, stBYM.model$waic$waic, stIntI.BYM.model$waic$waic))
)

row.names(dat.WAIC) = NULL

knitr::kable(dat.WAIC, caption = "WAIC of the fifferent models for dowry death cases") %>%  kable_styling(bootstrap_options = "striped", full_width = F, position = "center")

WAIC of the fifferent models for dowry death cases
model	WAIC
Spatial	659
SpatTemp no int	6757
SpatTemp typeI	6334

ANSWER: Although the spatial model has a lower WAIC it cannot be directly compared to the two spatio-temporal models, as it is essentially based on a different data set (i.e. data were aggregated removing the temporal component). The fair comparison in between the two spatio-temporal models, and we observe that the model with type I interaction performs better that the model without interaction.

Spatio-temporal analysis of dowry death crimes against women in Uttar Pradesh (India)

Hailey Gu

2023-03-29

1. Import and explore the data

2. Spatial model

3. Estimate the spatial fraction

4. Spatio-temporal model (no interaction)

5. Spatio-temporal model (type I interaction)

References